import matplotlib.pyplot as plt
import os,sys
import numpy as np
import pandas as pd
#evid1,evid2,o1,o2,tt1,tt2,dtcc,coeff,res_dt,net,station,pha
# 读取数据
df = pd.read_csv('./cc.csv')
print(f"数据形状: {df.shape}")
evid1=df['evid1']
evid2=df['evid2']

# 提取列数据
cc_dt = df['dt_cc']
pha = df['pha']
tt1 = df['tt1']
tt2 = df['tt2']
dt_cat = df['dt_cat']  # tt2-tt1
coeff = df['coeff']

# 创建相位掩码
maskarrayP = pha == "P"
maskarrayS = pha == "S"

print(f"P波数据量: {maskarrayP.sum()}, S波数据量: {maskarrayS.sum()}")

# 计算最小差值
min_diff = max(np.abs(cc_dt - dt_cat))
print(f"cc_dt - dt_cat 的最小值: {min_diff}")

# 修正后的循环
for cc_min in [0.5, 0.6, 0.8]:   
    print(f"\n=== 相关系数阈值: {cc_min} ===")
    
    # 方法1：分别处理P波和S波数据（推荐）
    # 提取P波数据
    df_P = df[maskarrayP].copy()
    df_S = df[maskarrayS].copy()
    
    # 应用相关系数筛选
    df_P_filtered = df_P[df_P['coeff'] > cc_min]
    df_S_filtered = df_S[df_S['coeff'] > cc_min]
    
    # 计算残差
    res_dt_P = df_P_filtered['dt_cc'] - df_P_filtered['dt_cat']
    res_dt_S = df_S_filtered['dt_cc'] - df_S_filtered['dt_cat']
 
    print(f"P波筛选后数据量: {len(df_P_filtered)}")
    print(f"S波筛选后数据量: {len(df_S_filtered)}")
    
    if len(res_dt_P) > 0:
        print(f"P波残差最小值: {max(np.abs(res_dt_P)):.6f}")
    else:
        print("P波无满足条件的数据")
        
    if len(res_dt_S) > 0:
        print(f"S波残差最小值: {max(np.abs(res_dt_S)):.6f}")
    else:
        print("S波无满足条件的数据")
        

    fig,axes=plt.subplots(1,2,figsize=(15,5))
    fig.suptitle(" with cc threshold > "+str(cc_min))
    ax1=axes[0]
    
    ax1.hist(res_dt_P,bins=50,edgecolor="black")
    text="mean {:.2} s \nstd= {:.2} s".format(np.average(res_dt_P),np.std(res_dt_P))
    ax1.set_title("P-wave")
    ax1.annotate(text,(0.6,0.8),xycoords="axes fraction")
    ax1.set_ylabel("Count")
    ax1.set_xlabel("CC_dt - Catalog_dt (s)")
    # ax1.plot([-15,15],[-15,15],"-r",linewidth=1)
    # ax1.axis("equal")
    # ax1.axis([-15,15,-15,15])
    # ax1.set_ylim(-15,15)
    # ax1.set_ylabel("Catalog_dt (s)")
    # ax1.set_xlabel("CC_dt (s)")
    ax2=axes[1]
    ax2.hist(res_dt_S,bins=50,edgecolor="black")
    ax2.set_title("S-wave")
    ax2.set_ylabel("Count")
    ax2.set_xlabel("CC_dt - Catalog_dt (s)")
    text="mean {:.2} s \nstd= {:.2} s".format(np.average(res_dt_S),np.std(res_dt_S))
    ax2.annotate(text,(0.6,0.8),xycoords="axes fraction")
    plt.tight_layout()
    fig.subplots_adjust(top=0.88)
    plt.savefig("./CC/ccthreshold_res_dt_{}.png".format(cc_min),dpi=300,inches_bbox="tight")
